# include "scene.h"
# include <stdio.h>
# include <stdlib.h>
# include <math.h>

# ifdef _OPENMP_
# include <omp.h>
# endif

extern FILE *file_open(FILE *lptr, const char *file, unsigned int t_cnt, unsigned int nfreq);
extern void conf_grid(struct grid_config *gc, struct grid *g, char *ecode);
extern void print_grid_conf(struct grid_config *gc, struct grid *g);
extern void grid_gen(struct grid_config *gc, struct grid *g, char *ecode);

extern short cmtrx_nthderv_CD(double *x, unsigned short nxpt, double *xderv, struct cmtrx *cm, double **tempmtrx[3]);

extern short cmtrx_1dim_interpol(double *x, unsigned short nxpt, double *xipl, unsigned short nxiplpt,\
 struct cmtrx *cm, double **tempmtrx);

// extern short cmtrx_multiply(struct cmtrx *mleft, struct cmtrx *mright, struct cmtrx *mres, double **tempmtrx);

static void cmtrx_distroy(struct cmtrx *cm);

# ifdef SOLVE_ENERGY
extern short heater_param(struct grid *g);
# endif

extern void mem_alloc(FILE *lptr, struct arr_struct *arr);
extern void mem_distroy(FILE *lptr, struct arr_struct *arr);

extern short solve(struct grid *g, struct arr_struct *depVar, struct arr_struct (*memFace)[3], struct arr_struct *memOut);

int main(int argc, char *argv[])
  {
   struct grid g;
   struct grid_config gc[Nsec];
   unsigned short isec, tn, nCPU, dim;
   short i, n;
   char ecode[3];
   struct arr_struct tnext, depVar[Nsec], memFace[Nsec][3], memOut;
   double *xmid;
   struct cmtrx *cm;
// short j;

   printf("Program starts: in 'main()'...\n");

#  ifdef LOG_FILE
   g.lptr = file_open(stdout, LOG_FILE, 0, 0);
   setlinebuf(g.lptr);
#  endif

#  ifdef TRACE_FILE
   g.trptr = file_open(stdout, TRACE_FILE, 0, 0);
   setlinebuf(g.trptr);
#  endif

   g.nTRD = nCPU = 1;

#  ifdef _OPENMP_
   fprintf(g.lptr, "Parallel computation using OpenMP.\n");

/* Setting number of threads more than the number of available physical CPUs. */
   g.nTRD = (nCPU = omp_get_num_procs()) + 0;
   omp_set_num_threads(g.nTRD);
/* Enabling dynamic and nested threading. */
   omp_set_nested(1);
   if(!omp_get_nested())
       {
        fprintf(g.lptr, "ERROR:: Nested threads are not enabled.\nExiting...\n");
        return -1;
       }

   fprintf(g.lptr, "Number of thread(s) used for this computation: %d\n", g.nTRD);
#  endif

   fprintf(g.lptr, "Number of available CPU(s): %d\n", nCPU);

/* --- Grid generation--- */
   conf_grid(gc, &g, ecode);
   if((ecode[0]|ecode[1]|ecode[2]) != 'o')
      {
       fprintf(g.lptr, "ERROR:: In 'main()': Failed to configure grid-parameters.\nError-code = %c:%c:%c.\nExiting...\n", ecode[0], ecode[1], ecode[2]);
       return -1;
      }

/* Printing grid-configuation parameters. */
   print_grid_conf(gc, &g);

/* Allocating arrays to store the grid-POINTS. */

/* Memory to store theta grid-points */
   g.th = (double *)malloc((g.nth + 2*NthG + 1)*sizeof(double));
   for(isec = 0; isec < Nsec; isec++)
      {
       g.r[isec] = (double *)malloc((g.nr[isec] + 2*NrG + 1)*sizeof(double));
       g.z[isec] = (double *)malloc((g.nz[isec] + 2*NzG + 1)*sizeof(double));
      }

   grid_gen(gc, &g, ecode);

   if((ecode[0]|ecode[1]|ecode[2]) != 'o')
      {
       fprintf(g.lptr, "ERROR:: In 'main()': Grid-generation failed.\nError-code = %c:%c:%c.\nExiting...\n", ecode[0], ecode[1], ecode[2]);
       return -1;
      }

/* --- Grid generation done --- */


/* ---- Determining co-efficient matrices for discretization ---- */
   fprintf(g.lptr, "Generating co-efficient matrices for dicretization and interpolation...\n");

/* Calculating the maximum number of grids. */
   n = g.nth+2*NthG;
   for(isec = 0; isec < Nsec; isec++)
       {
        if(g.nz[isec]+2*NzG > n) n = g.nz[isec]+2*NzG;
        if(g.nr[isec]+2*NrG > n) n = g.nr[isec]+2*NrG;
       }

/* Allocating temporary space to keep the mid-points. */
   xmid = (double *)malloc(sizeof(double)*n);

   fprintf(g.lptr, "\tin theta-direction...\n");

/* Determining discretization co-efficients for theta-direction. */
   for(i = 0; i < g.nth+2*NthG; i++) xmid[i] = 0.5*(g.th[i]+g.th[i+1]);

   cm = &g.dcmx_thcc2;
   cm->srcpos = 0.0;
   cm->respos = 0.0;
   cm->napprox = nDISCR;
   cm->nderv = 2;
   cm->srcpadi = 0;
   cm->respadi = 0;
   cm->srcpade = 0;
   cm->respade = 0;
   if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], g.nth+2*NthG, (double *)0, cm, (double ***)0) < 0) return -1;

   cm = &g.dcmx_thcc1;
   cm->srcpos = 0.0;
   cm->respos = 0.0;
   cm->napprox = nDISCR;
   cm->nderv = 1;
   cm->srcpadi = 0;
   cm->respadi = 0;
   cm->srcpade = 0;
   cm->respade = 0;
   if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], g.nth+2*NthG, (double *)0, cm, (double ***)0) < 0) return -1;

   cm = &g.dcmx_thccm;
   cm->srcpos = 0.0;
   cm->respos = 0.0;
   cm->napprox = 2;  //Yes, it's not a mistake. It'll be used to evaluate upwinding scheme.
   cm->nderv = nDISCR; //Yes, it's not a mistake. It'll be used to evaluate upwinding scheme.
   cm->srcpadi = 0;
   cm->respadi = 0;
   cm->srcpade = 0;
   cm->respade = 0;
   if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], g.nth+2*NthG, (double *)0, cm, (double ***)0) < 0) return -1;

   cm = &g.dcmx_thcs1;
   cm->srcpos = 0.0;
   cm->respos = 0.5;
   cm->napprox = nDISCR;
   cm->nderv = 1;
   cm->srcpadi = 0;
   cm->respadi = 1;
   cm->srcpade = 0;
   cm->respade = 1;
   if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], g.nth+2*NthG, &g.th[cm->respadi], cm, (double ***)0) < 0) return -1;

   cm = &g.dcmx_thsc1;
   cm->srcpos = 0.5;
   cm->respos = 0.0;
   cm->napprox = nDISCR;
   cm->nderv = 1;
   cm->srcpadi = 1;
   cm->respadi = 1;
   cm->srcpade = 0;
   cm->respade = 0;
   if(cmtrx_nthderv_CD(&g.th[cm->srcpadi], g.nth+2*NthG, &xmid[cm->respadi], cm, (double ***)0) < 0) return -1;

   cm = &g.dcmx_thss1;
   cm->srcpos = 0.5;
   cm->respos = 0.5;
   cm->napprox = nDISCR;
   cm->nderv = 1;
   cm->srcpadi = 1;
   cm->respadi = 1;
   cm->srcpade = 0;
   cm->respade = 0;
   if(cmtrx_nthderv_CD(&g.th[cm->srcpadi], g.nth+2*NthG, (double *)0, cm, (double ***)0) < 0) return -1;

   cm = &g.dcmx_thss2;
   cm->srcpos = 0.5;
   cm->respos = 0.5;
   cm->napprox = nDISCR;
   cm->nderv = 2;
   cm->srcpadi = 1;
   cm->respadi = 1;
   cm->srcpade = 0;
   cm->respade = 0;
   if(cmtrx_nthderv_CD(&g.th[cm->srcpadi], g.nth+2*NthG, (double *)0, cm, (double ***)0) < 0) return -1;

   cm = &g.icmx_thcs;
   cm->srcpos = 0.0;
   cm->respos = 0.5;
   cm->napprox = nDISCR-1;
   cm->nderv = -1; //Negative to show that this is not a derivative matrix.
   cm->srcpadi = 0;
   cm->respadi = 1;
   cm->srcpade = 0;
   cm->respade = 1;
   if(cmtrx_1dim_interpol(&xmid[cm->srcpadi], 2*NthG+g.nth, &g.th[cm->respadi], 2*NthG+g.nth-1, cm, (double **)0) < 0) return -1;

   cm = &g.icmx_thsc;
   cm->srcpos = 0.5;
   cm->respos = 0.0;
   cm->napprox = nDISCR-1;
   cm->nderv = -1; //Negative to show that this is not a derivative matrix.
   cm->srcpadi = 1;
   cm->respadi = 1;
   cm->srcpade = 0;
   cm->respade = 0;
   if(cmtrx_1dim_interpol(&g.th[cm->srcpadi], 2*NthG+g.nth, &xmid[cm->respadi], 2*NthG+g.nth-1, cm, (double **)0) < 0) return -1;

   fprintf(g.lptr, "\t...done.\n\tin r-direction...\n");

/* Determining co-efficient for discretization/interpolation in r-direction. */
   for(isec = 0; isec <= 1; isec++) // For sections 0 and 1.
     {
      fprintf(g.lptr, "\t\tfor isec = %d...\n", isec);
      for(i = 0; i < 2*NrG+g.nr[isec]; i++) xmid[i] = 0.5*(g.r[isec][i]+g.r[isec][i+1]);

      cm = &g.dcmx_rcc1[isec];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = 0;
      cm->respadi = 0;
      cm->srcpade = NrG-1;
      cm->respade = NrG-1;
      if( cmtrx_nthderv_CD(&xmid[cm->srcpadi], NrG+g.nr[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_rcc2[isec];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = nDISCR;
      cm->nderv = 2;
      cm->srcpadi = 0;
      cm->respadi = 0;
      cm->srcpade = NrG-1;
      cm->respade = NrG-1;
      if( cmtrx_nthderv_CD(&xmid[cm->srcpadi], NrG+g.nr[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_rccm[isec];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = 2; //Yes, it's not a mistake. It'll be used to evaluate upwinding scheme.
      cm->nderv = nDISCR; //Yes, it's not a mistake. It'll be used to evaluate upwinding scheme.
      cm->srcpadi = 0;
      cm->respadi = 0;
      cm->srcpade = NrG-1;
      cm->respade = NrG-1;
      if( cmtrx_nthderv_CD(&xmid[cm->srcpadi], NrG+g.nr[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_rcs1[isec];
      cm->srcpos = 0.0;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = 0;
      cm->respadi = 1;
      cm->srcpade = NrG-1;
      cm->respade = NrG;
      if( cmtrx_nthderv_CD(&xmid[cm->srcpadi], NrG+g.nr[isec]+1, &g.r[isec][cm->respadi], cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_rss1[isec];
      cm->srcpos = 0.5;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = 1;
      cm->respadi = 1;
      cm->srcpade = NrG-1;
      cm->respade = NrG-1;
      if( cmtrx_nthderv_CD(&g.r[isec][cm->srcpadi], NrG+g.nr[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_rss2[isec];
      cm->srcpos = 0.5;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 2;
      cm->srcpadi = 1;
      cm->respadi = 1;
      cm->srcpade = NrG-1;
      cm->respade = NrG-1;
      if( cmtrx_nthderv_CD(&g.r[isec][cm->srcpadi], NrG+g.nr[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_rsc1[isec];
      cm->srcpos = 0.5;
      cm->respos = 0.0;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = 1;
      cm->respadi = 1;
      cm->srcpade = NrG-1;
      cm->respade = NrG-1;
      if( cmtrx_nthderv_CD(&g.r[isec][cm->srcpadi], NrG+g.nr[isec]+1, &xmid[cm->respadi], cm, (double ***)0) < 0) return -1;

      cm = &g.icmx_rcs[isec];
      cm->srcpos = 0.0;
      cm->respos = 0.5;
      cm->napprox = nDISCR-1;
      cm->nderv = -1;
      cm->srcpadi = 0;
      cm->respadi = 1;
      cm->srcpade = NrG-1;
      cm->respade = NrG;
      if(cmtrx_1dim_interpol(&xmid[cm->srcpadi], NrG+g.nr[isec]+1, &g.r[isec][cm->respadi], NrG+g.nr[isec], cm, (double **)0) < 0) return -1;

      cm = &g.icmx_rsc[isec];
      cm->srcpos = 0.5;
      cm->respos = 0.0;
      cm->napprox = nDISCR-1;
      cm->nderv = -1;
      cm->srcpadi = 1;
      cm->respadi = 1;
      cm->srcpade = NrG-1;
      cm->respade = NrG-1;
      if(cmtrx_1dim_interpol(&g.r[isec][cm->srcpadi], NrG+g.nr[isec]+1, &xmid[cm->respadi], NrG+g.nr[isec], cm, (double **)0) < 0) return -1;

      fprintf(g.lptr, "\t\t...done.\n");
     }

//   For section 2.
      isec = 2;
      fprintf(g.lptr, "\t\tfor isec = %d...\n", isec);
      for(i = 0; i < 2*NrG+g.nr[isec]; i++) xmid[i] = 0.5*(g.r[isec][i]+g.r[isec][i+1]);

      cm = &g.dcmx_rcc1[isec];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = NrG-1;
      cm->respadi = NrG-1;
      cm->srcpade = NrG-1;
      cm->respade = NrG-1;
      if( cmtrx_nthderv_CD(&xmid[cm->srcpadi], 1+g.nr[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_rcc2[isec];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = nDISCR;
      cm->nderv = 2;
      cm->srcpadi = NrG-1;
      cm->respadi = NrG-1;
      cm->srcpade = NrG-1;
      cm->respade = NrG-1;
      if( cmtrx_nthderv_CD(&xmid[cm->srcpadi], 1+g.nr[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_rccm[isec];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = 2; //Yes, it's not a mistake. It'll be used to evaluate upwinding scheme.
      cm->nderv = nDISCR; //Yes, it's not a mistake. It'll be used to evaluate upwinding scheme.
      cm->srcpadi = 0;
      cm->respadi = 0;
      cm->srcpade = NrG-1;
      cm->respade = NrG-1;
      if( cmtrx_nthderv_CD(&xmid[cm->srcpadi], NrG+g.nr[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_rcs1[isec];
      cm->srcpos = 0.0;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = NrG-1;
      cm->respadi = NrG;
      cm->srcpade = NrG-1;
      cm->respade = NrG;
      if( cmtrx_nthderv_CD(&xmid[cm->srcpadi], 1+g.nr[isec]+1, &g.r[isec][cm->respadi], cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_rss1[isec];
      cm->srcpos = 0.5;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = NrG-1;
      cm->respadi = NrG-1;
      cm->srcpade = NrG-1;
      cm->respade = NrG-1;
      if( cmtrx_nthderv_CD(&g.r[isec][cm->srcpadi], g.nr[isec]+3, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_rss2[isec];
      cm->srcpos = 0.5;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 2;
      cm->srcpadi = NrG-1;
      cm->respadi = NrG-1;
      cm->srcpade = NrG-1;
      cm->respade = NrG-1;
      if( cmtrx_nthderv_CD(&g.r[isec][cm->srcpadi], g.nr[isec]+3, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_rsc1[isec];
      cm->srcpos = 0.5;
      cm->respos = 0.0;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = NrG-1;
      cm->respadi = NrG-1;
      cm->srcpade = NrG-1;
      cm->respade = NrG-1;
      if( cmtrx_nthderv_CD(&g.r[isec][cm->srcpadi], g.nr[isec]+3, &xmid[cm->respadi], cm, (double ***)0) < 0) return -1;

      cm = &g.icmx_rcs[isec];
      cm->srcpos = 0.0;
      cm->respos = 0.5;
      cm->napprox = nDISCR-1;
      cm->nderv = -1;
      cm->srcpadi = NrG-1;
      cm->respadi = NrG;
      cm->srcpade = NrG-1;
      cm->respade = NrG;
      if(cmtrx_1dim_interpol(&xmid[cm->srcpadi], 1+g.nr[isec]+1, &g.r[isec][cm->respadi], g.nr[isec]+1, cm, (double **)0) < 0) return -1;

      cm = &g.icmx_rsc[isec];
      cm->srcpos = 0.5;
      cm->respos = 0.0;
      cm->napprox = nDISCR-1;
      cm->nderv = -1;
      cm->srcpadi = NrG-1;
      cm->respadi = NrG-1;
      cm->srcpade = NrG-1;
      cm->respade = NrG-1;
      if(cmtrx_1dim_interpol(&g.r[isec][cm->srcpadi], g.nr[isec]+3, &xmid[cm->respadi], 1+g.nr[isec]+1, cm, (double **)0) < 0) return -1;

      fprintf(g.lptr, "\t\t...done.\n");

   fprintf(g.lptr, "\t...done.\n\tin z-direction...\n");


/* Determining discretization co-efficient in z-direction. */
   // For section-0 and 2;
      isec = 0;
CMTRX_ZSEC02:
      fprintf(g.lptr, "\t\tfor isec = %d...\n", isec);
      for(i = 0; i < g.nz[isec]+2*NzG; i++) xmid[i] = 0.5*(g.z[isec][i]+g.z[isec][i+1]);

      g.dcmx_zcc1[isec] = (struct cmtrx *)malloc(sizeof(struct cmtrx));
      cm = &g.dcmx_zcc1[isec][0];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = NzG-1;
      cm->respadi = NzG-1;
      cm->srcpade = 0;
      cm->respade = 0;
      if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], 1+g.nz[isec]+NzG, (double *)0, cm, (double ***)0) < 0) return -1;

      g.dcmx_zcc2[isec] = (struct cmtrx *)malloc(sizeof(struct cmtrx));
      cm = &g.dcmx_zcc2[isec][0];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = nDISCR;
      cm->nderv = 2;
      cm->srcpadi = NzG-1;
      cm->respadi = NzG-1;
      cm->srcpade = 0;
      cm->respade = 0;
      if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], 1+g.nz[isec]+NzG, (double *)0, cm, (double ***)0) < 0) return -1;

      g.dcmx_zccm[isec] = (struct cmtrx *)malloc(sizeof(struct cmtrx));
      cm = &g.dcmx_zccm[isec][0];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = 2; //Yes, it's not a mistake. It'll be used to evaluate upwinding scheme.
      cm->nderv = nDISCR; //Yes, it's not a mistake. It'll be used to evaluate upwinding scheme.
      cm->srcpadi = NzG-1;
      cm->respadi = NzG-1;
      cm->srcpade = 0;
      cm->respade = 0;
      if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], 1+g.nz[isec]+NzG, (double *)0, cm, (double ***)0) < 0) return -1;

      g.dcmx_zcs1[isec] = (struct cmtrx *)malloc(sizeof(struct cmtrx));
      cm = &g.dcmx_zcs1[isec][0];
      cm->srcpos = 0.0;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = NzG-1;
      cm->respadi = NzG;
      cm->srcpade = 0;
      cm->respade = 1;
      if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], 1+g.nz[isec]+NzG, &g.z[isec][cm->respadi], cm, (double ***)0) < 0) return -1;

      g.dcmx_zss1[isec] = (struct cmtrx *)malloc(sizeof(struct cmtrx));
      cm = &g.dcmx_zss1[isec][0];
      cm->srcpos = 0.5;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = NzG-1;
      cm->respadi = NzG-1;
      cm->srcpade = 0;
      cm->respade = 0;
      if(cmtrx_nthderv_CD(&g.z[isec][cm->srcpadi], 1+g.nz[isec]+1+NzG, (double *)0, cm, (double ***)0) < 0) return -1;

      g.dcmx_zss2[isec] = (struct cmtrx *)malloc(sizeof(struct cmtrx));
      cm = &g.dcmx_zss2[isec][0];
      cm->srcpos = 0.5;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 2;
      cm->srcpadi = NzG-1;
      cm->respadi = NzG-1;
      cm->srcpade = 0;
      cm->respade = 0;
      if(cmtrx_nthderv_CD(&g.z[isec][cm->srcpadi], 1+g.nz[isec]+1+NzG, (double *)0, cm, (double ***)0) < 0) return -1;

      g.icmx_zcs[isec] = (struct cmtrx *)malloc(sizeof(struct cmtrx));
      cm = &g.icmx_zsc[isec][0];
      cm->srcpos = 0.5;
      cm->respos = 0.0;
      cm->napprox = nDISCR-1;
      cm->nderv = -1;
      cm->srcpadi = NzG-1;
      cm->respadi = NzG-1;
      cm->srcpade = 0;
      cm->respade = 0;
      if(cmtrx_1dim_interpol(&g.z[isec][cm->srcpadi], 1+g.nz[isec]+1+NzG, &xmid[cm->respadi], g.nz[isec]+2, cm, (double **)0) < 0) return -1;

      g.icmx_zcs[isec] = (struct cmtrx *)malloc(sizeof(struct cmtrx));
      cm = &g.icmx_zcs[isec][0];
      cm->srcpos = 0.0;
      cm->respos = 0.5;
      cm->napprox = nDISCR-1;
      cm->nderv = -1;
      cm->srcpadi = NzG-1;
      cm->respadi = NzG;
      cm->srcpade = 0;
      cm->respade = 1;
      if(cmtrx_1dim_interpol(&xmid[cm->srcpadi], 1+g.nz[isec]+NzG, &g.z[isec][cm->respadi], g.nz[isec]+NzG, cm, (double **)0) < 0) return -1;

      fprintf(g.lptr, "\t\t...done.\n");

   if(isec == 0)
     {
      isec = 2;
      goto CMTRX_ZSEC02;
     }


   // For section-1.
      isec = 1;
//    for(isec = 0; isec < Nsec; isec++)
//      {
      fprintf(g.lptr, "\t\tfor isec = %d...\n", isec);
      for(i = 0; i < g.nz[isec]+2*NzG; i++) xmid[i] = 0.5*(g.z[isec][i]+g.z[isec][i+1]);

      g.dcmx_zcc1[isec] = (struct cmtrx *)malloc(2*sizeof(struct cmtrx));
      cm = &g.dcmx_zcc1[isec][0];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = 0;
      cm->respadi = 0;
      cm->srcpade = NzG-1;
      cm->respade = NzG-1;
      if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], NzG+g.nz[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_zcc1[isec][1];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = NzG-1;
      cm->respadi = NzG-1;
      cm->srcpade = NzG-1;
      cm->respade = NzG-1;
      if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], 1+g.nz[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;


      g.dcmx_zcc2[isec] = (struct cmtrx *)malloc(2*sizeof(struct cmtrx));
      cm = &g.dcmx_zcc2[isec][0];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = nDISCR;
      cm->nderv = 2;
      cm->srcpadi = 0;
      cm->respadi = 0;
      cm->srcpade = NzG-1;
      cm->respade = NzG-1;
      if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], NzG+g.nz[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_zcc2[isec][1];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = nDISCR;
      cm->nderv = 2;
      cm->srcpadi = NzG-1;
      cm->respadi = NzG-1;
      cm->srcpade = NzG-1;
      cm->respade = NzG-1;
      if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], 1+g.nz[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;


      g.dcmx_zccm[isec] = (struct cmtrx *)malloc(2*sizeof(struct cmtrx));
      cm = &g.dcmx_zccm[isec][0];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = 2; //Yes, it's not a mistake. It'll be used to evaluate upwinding scheme.
      cm->nderv = nDISCR; //Yes, it's not a mistake. It'll be used to evaluate upwinding scheme.
      cm->srcpadi = 0;
      cm->respadi = 0;
      cm->srcpade = NzG-1;
      cm->respade = NzG-1;
      if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], NzG+g.nz[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_zccm[isec][1];
      cm->srcpos = 0.0;
      cm->respos = 0.0;
      cm->napprox = 2; //Yes, it's not a mistake. It'll be used to evaluate upwinding scheme.
      cm->nderv = nDISCR; //Yes, it's not a mistake. It'll be used to evaluate upwinding scheme.
      cm->srcpadi = NzG-1;
      cm->respadi = NzG-1;
      cm->srcpade = NzG-1;
      cm->respade = NzG-1;
      if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], 1+g.nz[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;


      g.dcmx_zcs1[isec] = (struct cmtrx *)malloc(2*sizeof(struct cmtrx));
      cm = &g.dcmx_zcs1[isec][0];
      cm->srcpos = 0.0;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = 0;
      cm->respadi = 1;
      cm->srcpade = NzG-1;
      cm->respade = NzG;
      if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], NzG+g.nz[isec]+1, &g.z[isec][cm->respadi], cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_zcs1[isec][1];
      cm->srcpos = 0.0;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = NzG-1;
      cm->respadi = NzG;
      cm->srcpade = NzG-1;
      cm->respade = NzG;
      if(cmtrx_nthderv_CD(&xmid[cm->srcpadi], g.nz[isec]+2, &g.z[isec][cm->respadi], cm, (double ***)0) < 0) return -1;


      g.dcmx_zss1[isec] = (struct cmtrx *)malloc(2*sizeof(struct cmtrx));
      cm = &g.dcmx_zss1[isec][0];
      cm->srcpos = 0.5;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = 1;
      cm->respadi = 1;
      cm->srcpade = NzG-1;
      cm->respade = NzG-1;
      if(cmtrx_nthderv_CD(&g.z[isec][cm->srcpadi], NzG+g.nz[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_zss1[isec][1];
      cm->srcpos = 0.5;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 1;
      cm->srcpadi = NzG-1;
      cm->respadi = NzG-1;
      cm->srcpade = NzG-1;
      cm->respade = NzG-1;
      if(cmtrx_nthderv_CD(&g.z[isec][cm->srcpadi], g.nz[isec]+3, (double *)0, cm, (double ***)0) < 0) return -1;


      g.dcmx_zss2[isec] = (struct cmtrx *)malloc(2*sizeof(struct cmtrx));
      cm = &g.dcmx_zss2[isec][0];
      cm->srcpos = 0.5;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 2;
      cm->srcpadi = 1;
      cm->respadi = 1;
      cm->srcpade = NzG-1;
      cm->respade = NzG-1;
      if(cmtrx_nthderv_CD(&g.z[isec][cm->srcpadi], NzG+g.nz[isec]+1, (double *)0, cm, (double ***)0) < 0) return -1;

      cm = &g.dcmx_zss2[isec][1];
      cm->srcpos = 0.5;
      cm->respos = 0.5;
      cm->napprox = nDISCR;
      cm->nderv = 2;
      cm->srcpadi = NzG-1;
      cm->respadi = NzG-1;
      cm->srcpade = NzG-1;
      cm->respade = NzG-1;
      if(cmtrx_nthderv_CD(&g.z[isec][cm->srcpadi], g.nz[isec]+3, (double *)0, cm, (double ***)0) < 0) return -1;


      g.icmx_zsc[isec] = (struct cmtrx *)malloc(2*sizeof(struct cmtrx));
      cm = &g.icmx_zsc[isec][0];
      cm->srcpos = 0.5;
      cm->respos = 0.0;
      cm->napprox = nDISCR-1;
      cm->nderv = -1;
      cm->srcpadi = 1;
      cm->respadi = 1;
      cm->srcpade = NzG-1;
      cm->respade = NzG-1;
      if(cmtrx_1dim_interpol(&g.z[isec][cm->srcpadi], NzG+g.nz[isec]+1, &xmid[cm->respadi], NzG+g.nz[isec], cm, (double **)0) < 0) return -1;

      cm = &g.icmx_zsc[isec][1];
      cm->srcpos = 0.5;
      cm->respos = 0.0;
      cm->napprox = nDISCR-1;
      cm->nderv = -1;
      cm->srcpadi = NzG-1;
      cm->respadi = NzG-1;
      cm->srcpade = NzG-1;
      cm->respade = NzG-1;
      if(cmtrx_1dim_interpol(&g.z[isec][cm->srcpadi], g.nz[isec]+3, &xmid[cm->respadi], g.nz[isec]+2, cm, (double **)0) < 0) return -1;


      g.icmx_zcs[isec] = (struct cmtrx *)malloc(2*sizeof(struct cmtrx));
      cm = &g.icmx_zcs[isec][0];
      cm->srcpos = 0.0;
      cm->respos = 0.5;
      cm->napprox = nDISCR-1;
      cm->nderv = -1;
      cm->srcpadi = 0;
      cm->respadi = 1;
      cm->srcpade = NzG-1;
      cm->respade = NzG;
      if(cmtrx_1dim_interpol(&xmid[cm->srcpadi], NzG+g.nz[isec]+1, &g.z[isec][cm->respadi], NzG+g.nz[isec], cm, (double **)0) < 0) return -1;

      cm = &g.icmx_zcs[isec][1];
      cm->srcpos = 0.0;
      cm->respos = 0.5;
      cm->napprox = nDISCR-1;
      cm->nderv = -1;
      cm->srcpadi = NzG-1;
      cm->respadi = NzG;
      cm->srcpade = NzG-1;
      cm->respade = NzG;
      if(cmtrx_1dim_interpol(&xmid[cm->srcpadi], g.nz[isec]+2, &g.z[isec][cm->respadi], g.nz[isec]+1, cm, (double **)0) < 0) return -1;

      fprintf(g.lptr, "\t\t...done.\n");
//      }

   fprintf(g.lptr, "\t...done.\n");

   free(xmid);

   fprintf(g.lptr, "...done.\n");
/* -------- Generation of the co-effifient matrices are done. ------------ */


#  ifdef SOLVE_ENERGY
   fprintf(g.lptr, "NOTE: The energy equation will also be solved. Heater will be ON after %E seconds of the flow-start.\n", t_HEATER_ON);

/* Finding the heater's radius on the top surface. */
   if(heater_param(&g) < 0) return -1;
#  endif

/* Allocating memory. */
   for(isec = 0; isec < Nsec; isec++)
      {
    /* Allocating memory to store 3xp and 4x(u, v, w), and 1xT, while solving the Momentum and Energy equations. */
       depVar[isec].n1 = g.nr[isec] + 2*NrG;
       depVar[isec].n2 = g.nth + 2*NthG;
       depVar[isec].n3 = g.nz[isec] + 2*NzG;
       depVar[isec].nq = 15;
#      ifdef SOLVE_ENERGY
       depVar[isec].nq = 16;
#      endif
       mem_alloc(g.lptr, &depVar[isec]);

       g.u[isec] = depVar[isec].q[0];
       g.v[isec] = depVar[isec].q[1];
       g.w[isec] = depVar[isec].q[2];
       g.p[isec] = depVar[isec].q[3];
#      ifdef SOLVE_ENERGY
       g.T[isec] = depVar[isec].q[15];
#      endif

   /* Allocating memory to store A, B and C's at the cell-faces */
      for(dim = 0; dim < 3; dim++)
           {
            memFace[isec][dim].n1 = g.nr[isec] + (3-dim)/3;
            memFace[isec][dim].n2 = g.nth + dim%2;
            memFace[isec][dim].n3 = g.nz[isec] + (1+dim)/3;
            memFace[isec][dim].nq = nRK;
            mem_alloc(g.lptr, &memFace[isec][dim]);
           }
      }

/* Allocating memory to store u,v,w,T at the outlet for tn_BACKUP number of time-steps.
   It is required for application of the Orlansky BC. */
   memOut.n1 = g.nr_gap;
   memOut.n2 = g.nth;
   memOut.n3 = zCells_BACKUP;
   memOut.nq = (DEPVAR_N-1)*tn_BACKUP;
   mem_alloc(g.lptr, &memOut);

   for(tn = 0; tn < tn_BACKUP; tn++)
     {
      g.outBC.u[tn] = memOut.q[tn];
      g.outBC.v[tn] = memOut.q[tn_BACKUP+tn];
      g.outBC.w[tn] = memOut.q[2*tn_BACKUP+tn];
#     ifdef SOLVE_ENERGY
      g.outBC.T[tn] = memOut.q[3*tn_BACKUP+tn];
#     endif
     }

/* Allocating memory to store dependent variables at outlet-cell. It will be used in application of the Orlansky BC. */
   tnext.n1 = DEPVAR_N-1; /* Corresponding to u, v, w, T */
   tnext.n2 = g.nr_gap;
   tnext.n3 = g.nth;
   tnext.nq = 1;
   mem_alloc(g.lptr, &tnext);
   g.outBC.tnext_ug = tnext.q[0][0];
   g.outBC.tnext_vg = tnext.q[0][1];
   g.outBC.tnext_wg = tnext.q[0][2];
#  ifdef SOLVE_ENERGY
   g.outBC.tnext_Tg = tnext.q[0][3];
#  endif

/* Memory allocation done. */

/* Solving the NS-Equation. */
   solve(&g, depVar, memFace, &memOut);

/* Freeing memory. */
   for(isec = 0; isec < Nsec; isec++)
      {
       cmtrx_distroy(&g.dcmx_rcc1[isec]);
       cmtrx_distroy(&g.dcmx_rcc2[isec]);
       cmtrx_distroy(&g.dcmx_rcs1[isec]);
       cmtrx_distroy(&g.dcmx_rss1[isec]);
       cmtrx_distroy(&g.dcmx_rss2[isec]);

       mem_distroy(g.lptr, &depVar[isec]);
       for(dim = 0; dim < 3; dim++) mem_distroy(g.lptr, &memFace[dim][isec]);
       free(g.r[isec]);
       free(g.z[isec]);
      }
   cmtrx_distroy(&g.dcmx_thcc2);
   cmtrx_distroy(&g.dcmx_thcs1);
   free(g.th);
   mem_distroy(g.lptr, &memOut);
   mem_distroy(g.lptr, &tnext);


#  ifdef TRACE_FILE
   if(g.trptr != stderr) fclose(g.trptr);
#  endif

   fprintf(g.lptr, "...'main()' ends!\n");
   if(g.lptr != stderr) fclose(g.lptr);

   printf("...program ends.\n");

   return 0;
  }

static void cmtrx_distroy(struct cmtrx *cm)
   {
    unsigned short i;

    for(i = 0; i < cm->nrows; i++) free(cm->mtrx[i]);
    free(cm->mtrx);
    free(cm->istart);
    free(cm->iend);

    return;
   }
